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Abstract 


Cross stream modes are defined as solutions to the Orr-Sommerfeld equa- 
tion which are propagating normal to the flow direction. These modes are 
utilized as a basis for a Hilbert space to approximate the spectrum of the Orr- 
Sommerfeld equation with plane Poiseuille flow. The cross-stream basis leads 
to a standard eigenvalue problem for the frequencies of Poiseuille flow insta- 
bility waves. The coefficient matrix in the eigenvalue problem is shown to 
be the sum of a real matrix and a negative-imaginary diagonal matrix which 
represents the frequencies of the cross stream modes. The real coefficient 
matrix is shown to approach a Toeplitz matrix when the row and column 
indices are large. The Toeplitz matrix is diagonally dominant, and the di- 
agonal elements vary inversely in magnitude with diagonal position. The 
Poiseuille flow eigenvalues are shown to lie within GerSgorin disks with radii 
bounded by the product of the average flow speed and the axial wavenumber. 
It is shown that the eigenvalues approach the GerSgorin disk centers when 
the mode index is large, so that the method may be used to compute spectra 
with an essentially unlimited number of elements. When the mode index is 
large, the real part of the eigenvalue is the product of the axial wavenumber 
and the average flow speed, and the imaginary part of the eigenvalue is iden- 
tical to the corresponding cross stream mode frequency. The cross stream 
method is numerically well-conditioned in comparison to Chebyshev based 
methods, providing equivalent accuracy for small mode indices and superior 
accuracy for large indices. 



1 Introduction 


The Orr-Sommerfeld equation represents the small perturbations of other- 
wise parallel incompressible flow. As such, it has received a great deal of 
study. Most of these results are summarized in the excellent texts by Chan- 
drasekhar [1] and by Drazin and Reid [2]. The primary objective of most 
investigations has been to define stability boundaries. These boundaries can 
be defined either by study of the growth and decay of modes in time or by 
study of their growth and decay in space. In either approach, it has been 
found advantageous to cast the problem in the form of a linear eigenvalue 
problem. Chandrasekhar [1], Dolph and Lewis [3], and Grosch and SaJwen [4] 
used transcendental basis functions from fourth order operators to estimate 
eigenvalues of the Orr-Soinmerfeld equation. Orzag [5] utilized Chebyshev 
polynomials in a study of the temporal eigenvalue problem, while Bridges and 
Morris [6] applied these polynomials in a study of the spatial formulation. 

The movement toward Chebyshev polynomials was motivated by their at- 
tractive mathematical properties and by the disappointing results obtained 
from earlier approaches. Dolph and Lewis [3] had computed temporal eigen- 
values using a set of basis functions derived from the Orr-Sommerfeld equa- 
tion itself. Their rationale was that these functions would provide a natural 
basis for the problem at hand and could thus be expected to provide good 
approximations to the general Orr-Sommerfeld equation. The basis utilized 
by Dolph and Lewis was the countable set of solutions to the Orr-Sommerfeld 
equation for modes propagating perpendicular to the steady flow direction. 
For these modes, the Orr-Sommerfeld equation becomes a fourth-order or- 
dinary differential equation with constant coefficients and thus has exact 
solutions given by the elementary transcendental functions. 

But it is clear that Dolph and Lewis [3j were not entirely satisfied with 
their results. This may have been, in part, because they were still struggling 
with the numerics of the eigenvalue solution process, which is almost taken 
for granted today. The purpose of this paper is to re-examine the application 
of the basis functions of Dolph and Lewis [3] to the problem of determining 
the temporal stability characteristics of the Orr-Sommerfeld equation. Plane 
Poiseuille flow will be used for this study. 

The characteristic solutions of the Orr-Sommerfeld equation are defined 
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by the system 


L 4 (f) — uL 2 (j> 

«±0 = o 

*'(± 1 ) = 0 

where L 4 is a fourth order operator and L 2 is a second order operator. Let 
{^ 0 ) V'u- • ■} be a basis for a Hilbert space with inner product < u,v >. A 
solution is defined by the expansion 


OO 

<t>=Yl 

n=0 

which gives the matrix eigenvalue problem 

(ipm, Ulpn) a n = U (ipm, L 2 Tp n ) a n 

whose solutions give eigenvalues and associated eigenvectors a n ». This 
paper compares two methods of solving the above eigenvalue problem for the 
case of Poiseuille flow. The first method utilizes the Chebyshev basis with 
the Lanczos tau criterion, and the second utilizes the cross-stream basis of 
Dolph and Lewis [3] with the Galerkin criterion. 

Section 2 restates the Orr-Sommerfeld equation to define fully the opera- 
tors and other notation. The Chebyshev/Tau method is applied to Poiseuille 
flow in section 3 and it is demonstrated that this method suffers from poor 
numerical conditioning. Section 4 develops the cross stream basis in detail 
and gives the procedures for evaluation the basis functions to machine pre- 
cision. The cross stream formulation of the eigenvalue problem is given for 
general flows in section 5 and applied to the case of Poiseuille flow in section 
6. The summary and conclusions are given in section 7. 


2 


2 The Orr-Sommerfeld Equation 

The stability of incompressible flow in a channel is defined by characteristic 
solutions of the Orr-Sommerfeld equation. This paper presents an analysis 
of that equation emphasizing the special case of Poiseuile flow where the 
density and viscosity are constant. Figure 1 shows the flow in a channel. We 
use the same coordinates as in Drazin and Reid [2] since this is the most 
comprehensive and current work on the subject. In this system, x is the 
flow direction, y is horizontal and perpendicular to the flow, and z is vertical 
and perpendicular to the flow. The origin of coordinates is placed midway 
between the plates so that the lower and upper walls are located at z = ±h. 



Poiseuille flow is formed by having an axial pressure gradient to induce 
the flow. The x-momentum equation for the steady flow gives 


dp dr xz d ( dU\ 


0 ) 


where p is the viscosity. Since the shear stress gradient t' xx is constant, the 
stress varies linearly across the channel. The steady flow is then a parabola 


3 


( 2 ) 


and the maximum velocity is in the center of the channel. 

U(z) = U max (l- ^ 

The Orr- Sommer feld equation for velocity perturbations in the vertical di- 
rection (Drazin and Reid, [2, p. 156]) for a general velocity profile U(z ) 
is 

(U r - Cr)(w" - k^w) + (w"" - 2k^w" + k*w) - U"w = 0 (3) 

ICr 

where v = p/p is the kinematic viscosity. 

The symbols in equation 3 are explained below in figure 2. The figure 



Figure 2: Coordinates and wavenumber conventions. 


shows the wave vector geometry, which is taken according to a standard 
spherical coordinate convention. 

k = fcjCj "l - kyGy "I” kz&z 

kr = ks'mO 
k z = k cos 0 

k - Jk* + (4) 
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k x 

ky 

kr 


k r cos <j) 


hr sin 0 

•fii + 


Here K is the magnitude of the radial (horizontal) wave vector, which makes 
an angle <f> with respect to the downstream flow velocity, that is 


<f> = arctan 



The polar angle 0 is defined by the components k r and k z . 


(5) 


9 = arctan — (6) 

K z 

Here we have departed from Drazin and Reid is using k x and ky in place of 
their wavenumber symbols a and j3. We also depart from their convention in 
showing the explicit dependence of the variables on the propagation angle <f>. 
The component of U in the direction of the radial wave vector is U r = U cos <f>. 

The vertical component of the perturbation velocity is assumed to be a 
wave of the form 

We^r ( C08 ^ n <f> y—Crt) ^ 

If we define a radial position vector 


r = e x x + eyjj 


( 8 ) 


then the fluctuating vertical velocity can be given in compact form 

we ifc?-u>t) ( 9 ) 

Here, Cr is the complex radial propagation speed of the wave. The frequency 
of the wave is therefore 

u — k r Cr — k r cos <j>c— k x c (10) 

where c = Cx is the complex downstream speed used by Drazin and Reid. 
We use nondimensional variables based on the maximum velocity Umax , and 
the channel half-width h. The non-dimensional frequency is then u/hfUmax, 
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the non-dimensional wave number is k r h, and the nondimensional kinematic 
viscosity is ^/(Umaxh), which is the same as the inverse Reynolds number 
R~ l . The form of equation 3 remains the same, but we rewrite it as 

ivL A w + krCos<f>(u L?w - U"w ) = ujL 2 w (11) 

where the operators are defined as follows 

‘’-(s?-*) < 12 > 

t4= (^- 2 ^ + 4 < 13 > 

For Poiseuille flow, the nondimensional flow velocity function and its second 
derivative would be 


U={ 1 - z 2 ) 
U" = -2 


(14) 

(15) 


3 Chebyshev/Tau Spectral Method 

Orzag utilized a basis of Chebyshev polynomials with the Lanczos tau cri- 
terion to obtain accurate spectral solutions to the Orr-Sommerfeld equation 
for plane Poiseuille flow. This method is reviewed quickly here and Orzag’s 
results are replicated as a basis for comparison to the cross-stream spectral 
method. 

3.1 Review of the Chebyshev/Tau Method 

The Chebyshev method utilizes a series representation of w(z) in terms of 
the Chebyshev polynomials of the first kind. 

OO 

w(z) - ^2 a n T n (z) (16) 

n=0 

It was shown by Orzag that the derivatives of w(z) can be given by series of 
the same form as equation 16, provided that the coefficients in the series for 
the derivatives are defined in terms of the coefficients of the velocity function 
w. That is, the second derivative is given by the series 

OO 

w"(z) = £ afT n (z) (17) 

n=0 

where the coefficients in the series for the second derivative are defined in 
terms of the coefficients in the series for the function, i.e. 

1 00 

<*£? = — n (™ 2 “ ™ 2 )on, m> 0 (18) 

^Tl n— m+2 

n=m(niod 2) 

where 

(0, if m < 0; 

Cm = | 2, if m = 0; (19) 

1 1, if m > 0. 

Similar formulas for the fourth derivative are 

w""{z) = fx 4 >r n (z) (20) 

n=0 
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o£ } = ob~ £ P(m,n)a n> m> 0 (21) 

^ C tn n=m|4 

n~m (mod 3) 

P(m, n) = n [n 2 (n 2 - 4) 2 - 3m 2 n 2 (n 2 - m 2 ) - m 2 (m 2 - 4) 2 ] (22) 

These formulas are easily coded, but contain the origin of a significant prob- 
lem of numerical conditioning. When placed in a matrix format by truncating 
the series to a finite number of terms, say a polynomial of degree M - 1-1, the 
largest element in the matrix representing the fourth derivative operator is of 
order M 7 . The derivative operator matrices are upper triangular, with zeros 
on the diagonal, and the matrix equivalent of the operator L 4 will have di- 
agonal elements fc 4 , so that the smallest non-zero elements (on the diagonal) 
are of order unity, while the largest (in the upper right corner) are of order 
M 7 . 

In the case of Poiseuille flow, the only other relation needed is the ex- 
pression for the function z 2 w(z). The formula for this term is found from 
known recurrence formulas for Chebyshev polynomials. As given by Orzag, 
the relation is 

1 00 

Z w(z) = — ) ’ [Cn-2 a n-2 + ( Cn 4" C, j— l) &n 4" U«+ 2 ] (23) 

q n=0 

The matrix equivalent of this series is a tri-diagonal matrix operator with 
simple fractional elements such 1/2 and 1/4. This matrix is designated here 
by the symbol [Z 2 \. 

Using the symbols [L 2 ] and [L 4 ] to designate the second and fourth order 
operators, respectively, the matrix equivalent of the Orr-Sommerfeld equa- 
tion for Poiseuille flow is 

k, [21/] + [[/] - [Z 2 ]] (£ 2 |] { a„ } + iv[L 4 \ { } = u,[i, 2 ] { a „ } (24) 

In the Lanczos tau method, the rows of the above matrix equation rep- 
resenting the highest-degree polynomials are replaced with the constraints 
representing the boundary conditions on w and it/. In this way, the bound- 
ary conditions are satisfied exactly by the solutions to the matrix eigenvalue 
problem. 

Equation 24 has the form of a generalized eigenvalue problem, that is 

Ax = \Bx (25) 
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Conversion of this equation to a standard eigenvalue problem would require 
a matrix inverse. While this is generally possible for matrices of reasonable 
size, the choice is usually to use direct methods, such as the QZ algorithm, for 
a direct solution of the generalized eigenvalue problem. The results presented 
below are from the IMSL [8] library routine GVCCG for the complex gen- 
eralized eigenvalue problem. Measures of the accuracy of the computations 
were generated by the IMSL routine CPICG. 

3.2 Results of the Chebyshev/Tau Method 

Computations were made using the Chebyshev basis to find eigenvalues of 
the first 32 even modes, that is, where n = {0, 2, ■ • • 62}. Results are shown 
in figure 3 for two computational word sizes, the 32-bit word and the 64- 
bit word. Points shown in the figure have been rotated through 90° in the 
complex plane by plotting iu instead of u. The first quadrant then con- 
tains eigenvalues of stable modes. The word size clearly has a strong effect 



Figure 3: Computed eigenvalues iu from the Chebyshev/Tau method. K = 
1.0, (j> = 0.0, Re = 10 4 , and 32nd order matrix. Gray symbols show results 
of 32-bit words, and black show results of 64-bit words. 

on the results of the computation. The eigenvalues predicted by the 32-bit 
computation bear little or no relation to those predicted by the 64-bit com- 
putation, so that the less accurate computations must be incorrect. Accurate 
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computations with the Chebyshev basis require a large computer word. This 
is because the Chebyshev formulation is numerically ill-conditioned, that is, 
the condition numbers of the matrices in the generalized eigenvalue prob- 
lem are large. The IMSL subroutine CPICG issued warnings that the 32-bit 
computations shown in figure 3 would be inaccurate. 

There is another problem apparent in figure 3. One normally expects 32 
eigenvalues from a 32nd order matrix, and a careful count of the points in 
figure 3 shows only 25 black symbols and 23 gray. Unless there are repeti- 
tions, the others must be off-scale, but where? The answer is that they are 
scattered hither and yon in the complex plane, far from the origin, and in 
all four quadrants. The black symbol just to the left of the imaginary axis, 
near the point = 2.2, is normally called the“critical” eigenvalue, or 

fastest growing mode, but an examination of all the eigenvalues produced 
by this computation would indicate modes which grow much faster. Similar 
results are found also for the subcritical Reynolds numbers Re < 5772. It 
would appear that researchers who report results of “critical eigenvalue” and 
“critical Reynolds number” computations from the Chebyshev/Tau method 
are exercising judgements about the spectra in their conclusions. While 
these opinions are no-doubt carefully considered and responsible, the need 
for such judgements presents a danger which becomes more significant when 
computations are made for problems which have been explored to a lesser 
extent. It is highly desirable to have a spectral method which is relatively 
insensitive to computer-word size and which produces spectra such that the 
least-attenuated mode can be defined absolutely by the maximum imaginary 
part of the eigenvalue. This paper will show that the cross-mode method 
achieves these computational goals. 
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4 Cross-Stream Modes 


Dolph and Lewis [3] used solutions of the Orr-Sommerfeld equation with 
<l> = ±?r/2 as a basis for a spectral approximation. In this case, oos<f> = 0, 
and the Orr-Sommerfeld reduces to 

iuL^w = uL 2 w (26) 

The physical interpretation of these solutions is that the wave vector k r is 
pointing across the flow direction, so that these solutions will be called cross- 
stream modes here. 

4.1 Even Modes 

An even solution is 


t \ a 

w(z) = - 

where k is the spherical wavenumber and a is a constant to be defined later. 
This solution clearly satisfies the boundary conditions iu(±l) = 0. The 
boundary conditions it/ (±1) = 0 are then satisfied if 

0 > k z tan k z = -kr tanh k r > -oo, 0 < k r < oo (28) 

The function k z tan k z in equation 28 is even and, when k z is small, ap- 
proaches k 2 . Since kr tanh k T is even also, it approaches k 2 when k r is small. 
But there is no “small-/:*” solution to equation 28 because k z tan k z is posi- 
tive while — k r tanh kr is negative. The first solution k z o > 0 is thus found in 
the interval (ir/2 < k z o < n). Since k z tan k z is monotonically increasing in 
this interval, there is one and only one solution within the interval. Since the 
solutions k zn to equation 28 are real positive numbers, the following change 
of variable is introduced to facilitate their computation. 


cos k z z - cos k : 


cosh krZ 
coshAv 


(27) 


kin^kr) — K n 6 n (k r ) (29) 
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Then equation 28 becomes 

tan<5„ = — 3-7-, /?n(Av) = AvtanhAv, n = 0,2,4,- •• (31) 

This equation has a unique solution 0 < 6 n (k r ) < ir/2 for each selected 
even index n. Thus, equation 28 has a countable set of solutions which 
may be labled by the set of non-negative even integers. When Av is large, <5 n 
approaches 7r/2, and the solutions k zn approach (n+l)7r/2 from above. When 
kr is small, S n (k r ) approaches 0, and the solutions k xn approach (n + 2)tt/2 
from below. 

Asymptotic formulas for the solutions 8 n {K ) of equation 31 can be found 
by writing it as an arctangent function of the parameter Kn- 


S(f n ) = arctan 




A Taylor series in the parameter e n then gives the desired formula. 


0n (3-/Uff (30 - 20/7 n + 

«n 15«® 


(32) 

(33) 


Although it is valid in concept only for large n, when Kn » fc r > this formula 


is valid to four digits even in the case n = 0, Av = 1.0, kq = n. It thus 
provides a good initial estimate for a numerical solution of equation 32 for 
moderate values of k r . 

Numerical solutions for the wavenumbers are based on equation 32, writ- 
ten as a function of a variable 0 < £ = 28 n /ir < 1 together with its first 


derivative. 


where 



Figure 4 shows these functions for n = 0, 1 , 2, 00. The functions are nearly 
straight lines, so that the zero-crossing points £o(ri) are easily found numer- 
ically. The slopes of the curves are positive and nearly unity, so that the 
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errors in the functions axe roughly the same as the error in the solution. 
The solutions are found by starting with the estimate of the asymptotic for- 
mula 33 and then iterating with Newtons method using equations 34, 35 
until a desired accuracy is reached. Using a 64-bit computer word, which 
is equivalent to about 15 digits accuracy, the numerical solutions have been 
computed with an error no greater than 10 -14 . 


/n(0 



0 0.25 0.5 0.75 1 

£ 

Figure 4: Graphical solution for vertical wavenumbers, 7 „ = 1. 


4.2 Odd Modes 

An odd solution to equation 26 is 


The characteristic equation is again given by the boundary conditions on vJ 
1 > t” tan k z — -j- tanh fc r > 0, 0 < Ay < oo (38) 

Kz rCr 

Equation 38 also has a countable set of solutions. Since (tan k z )/{k z ) > 1 
in the interval ( — tt/2, 7t/ 2) and since (tanhAy)/Ay < 1 in the interval 0 < 
kr < oo, there is no solution within the central interval (— tt/2 < k z < 7t/2). 


sin k z z — sin k. 


sinh k r z 
sinh k r 


(37) 
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However, there is a single solution within each interval ((n 4- 1 )tt/2 < k zn < 
(n + 2)7r/2, n = 1, 3,5 •• -.Therefore, equation 38 has a countable set of 
solutions which may be labeled by the positive odd integers. When K is 
large, these solutions approach (n 4- 1)7t/2 from above (but now n is odd), 
and, when fcr is small, they approach a limit from below, but the limit is 
less than (n 4- 2)7r/2. Using the same ofTset variable definition as in the case 
of even modes, but with odd indices n, equation 38 is written in arctangent 
form. 

6(e n ) = arctan 

This equation is identical to the one for the even mode wavenumbers, except 
for the definition of /3 n , so that the asymptotic formula 33 and the function 
for numerical solution, equation 34 are applicable to the odd modes as well 
as the even. 


Priori 

1 - e n <5(e n ) 


0 n (k r ) = kr cothAv, n=l,3,5, (39) 


0 



0 8 16 24 32 40 48 56 64 

n 

Figure 5: Wavenumber offsets for 64 modes, kr = 1. 

Figure 5 shows the wavenumber offsets <5„ as computed by the asymp- 
totic formula. The points are plotted as log 10 6 n versus n. There are two 
patterns — the “upper” points are the odd indices, and the ‘lower” are the 
even. This is because the function 0 n (kr) is always larger for the odd indices. 
These offset variables become small quickly as the mode index increases. 
Consequently, the small-angle approximations for sinfi„ and cos 6 n will be 
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accurate for all but possibly the “lowest” modes. Figure 6 shows the accu- 
racy of the asymptotic formula for the offsets, compared to the numerical 
solution of the offset equation by Newton’s method. The asymptotic formula 



n 


Figure 6: Accuracy of asymptotic formula for wavenumber offsets for 64 
modes, k T — 1. 

has an absolute error less than 10“ 12 for mode indices greater than about 48. 
The formula is more accurate for the even indices than for the odd, but is an 
excellent approximation in either case. 

4.3 Combined Modes 

4.3.1 Definitions and identities 

The combined set of all real-valued solutions to equation 28 and to equa- 
tion 38 is thus a countable set of distinct wavenumbers which increase mono- 
tonically with the index n. 

. 7T - . 3n _ . 

{ 2 < kzo < 7T < k zl < — <k z2 < 27T < • • •} (40) 

Associated with these wavenumbers are characteristic solutions, called cross- 
modes here. It will be helpful to utilize the following identities in the defini- 
tion of these modes. 

cos k zm = — (— l)t m / 2 l cos 6 m , m even (41) 
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m odd 


sin k zm = -(-l) Im/21 cos<5 m , 


(42) 


The cross-modes are then defined as linear combinations of two “mode parts” 
) and i> n (z). 


w n (z) = 

^ [<f>n(z) + (-l) [n/2] cos M»„(z)] n = 0, 1, 2 • • • 

(43) 

Mz) = 

f cos k zn z, n even 
\ sin k zn z, n odd 

(44) 

i>n(z) = 

( cos ikrz/ oos ikr, n even 
1 sin ik r zj sin ikr, n odd. 

(45) 


The function ip n (z) is usually expressed in terms of real hyperbolic functions, 
but it will be convenient later to use the complex forms shown in equations 45. 
Note that the division operation [n/2] in equation 43 is an integer operation 
with an integer result. 

It is easy to show that 


= o 

(46) 

= -klMz) 

(47) 

L w n (z) — — Qf n fc n <^ n (z) 

(48) 


These results will be useful later in the development of spectral approxima- 
tions to Orr-Sommerfeld equation. 


4.3.2 Inner products, norms, and orthogonality 
An inner product is defined as an integral over the interval [-1,1]. 

if, 9) = \j_ i f(z)9(z)dz (49) 

The natural norm of a function is the square root of the inner product of the 
function with itself. 

ll/ll = </, /) 1/2 (50) 

The constants a n can be defined using the above definitions and the or- 
thonormal condition of Dolph and Lewis [3] 

^ j [k?w m (z)w n (z) + w' m {z)w' n {z)} dz = 6 mn (51) 
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Integrating equation 51 by parts and using the definition 49 gives the identity 

{^m j L Wn) — S m n (52) 

The cross-stream modes w n (z ) are not orthogonal; however, they are “partly 
orthogonal”, meaning that their two parts, as given by equation 41, are 
orthogonal. More precisely, the orthogonality properties of the two parts are 

(<Pm,i>n) = 0, all m,n (53) 

Since the functions ip are either even or odd, depending on the parity of the 
index, there is also the condition 


{ipm,ipn) = o, m + n odd (54) 

It can also be shown, using the characteristic equations for k zn , that the 
functions <p n are orthogonal. 


&rnn (4*n i 0n) 


(55) 


Combining this result with equation 52 defines the constants a* through the 
following equation. 

= (<f>n,(f > n )~ l/2 = ll^nll -1 (56) 

The inner product {(p n ,<p n ) is given by 

(<pn, (pn) = ^ (1 + (- l) n sine 2 k xn ) (57) 

where the “sine” function is defined as 


sin a: 

sine x = 

x 

Similarly, the inner product of the function ip n with itself is 


with 


(V'n.V'n) 


1 + (— l) n sinch2Av 
1 + (— l) n cosh2Av 


. , sinh x 

sinch x = 

x 


(58) 


(59) 

(60) 
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( 61 ) 


Using the expression for the vertical wavenumbers gives 



and the norm of the function 4> n is 



Since <5 n — » 0 for large n, it is clear that the above inner product approaches 
1/2. The final formula for the constants is then 



(63) 


Note that a n = 0(V 2) for all values of kr and n. The norm ||0„|| is a positive 
real number, roughly y/2/2 for any k r> including the limits where k T — * 0 and 
kr — ► oo. The norm of the function ip n is 


IUII- 1 /sinh2fc r ±2fc r V /2 / + ifn 
M ^ n " \ cosh2A: r ±l ) ’ l- ifn 


if n is even, 
is odd. 


(64) 


The norm ||^„|| approaches a positive limit, 1 for even modes and \/3/3 for 
odd, as kr approaches zero and approaches zero in the limit of large kr. The 
cross-mode functions are defined finally using the norm of the functions <f> n . 


4 >n + (~ l) |n/2] COS 6 n 1 p n 

K ||0n|| 


The norm of the cross-mode w n is 


lk.ll 


(ll^n|| 2 + COS 2 6n||^ l nH 2 ) 1 ^ 2 
^n||<Ani| 


(65) 


( 66 ) 


These results show that the cross- modes, as defined here, approach 0 uni- 
formly in z when either kr or n becomes large, because the magnitude of 
the spherical wavenumber k n becomes large. The derivatives of the cross- 
modes are finite, however, for these cases. The Dolph-Lewis normalization 
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condition, equation 51 can be given in terms of the norms ||ti>„|| and ||u£|| 
as 

IKII* + *? IKII 2 = 1 (67) 

This is the equation for an ellipse with ||iw n || plotted as the abscissa and 
IKII plotted as the ordinate. The “semi-major” axis is A:” 1 . , and the 
“semi-minor” axis is 1. Of course, if Av > 1 the major axis becomes the 
ordinate instead of the abscissa. But this equation makes clear the following 
bounds on the norms of w n and w' n . 


IKII < k r 1 (68) 

IKII < 1 (69) 

There is also the following explicit equation for the norm of the derivative in 
terms of the norms of the mode parts. 



*£Jl0n|| 2 - A:?cos 2 6 n l|V> w || 2 

KU nil 2 


(70) 


The ratio k tn /k n approaches unity for large k zn , consequently, for fixed Av 
and large index n, the norms have the following asymptotic limits. 


IKII ~ 0, n — ► oo (71) 

IKII ~ 1, n-» oo (72) 

Figures 7 and 8 illustrate the cross-modes and the norm properties de- 
scribed above. Figure 7 is the lowest mode, n = 0. Only the domain 
[0 < z < 1] is shown because the mode is even. The solid curve is the 
mode function w n (z) and the dashed curve is its derivative w' n {z). Note that 
both curves approach zero at the boundary as required by the boundary con- 
ditions. The norms of the function and its derivative are similar in magnitude 
for this mode. 

Figure 8 shows the eighth cross-mode. It is clear that the mode has 
smaller norm than its derivative norm in this case. The derivative of the 
mode is roughly a sine function. The norm of the derivative, which is the 
rms value taken over the channel width, is near unity for this mode. 

The cross-modes are non-orthogonal, or non-normal. A measure of rel- 
ative normality of different modes is the “angle” ') rnn between two different 
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Figure 8: Cross mode Wg(z) and its derivative u4(*)> K = 1. 
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modes. This angle can be defined as 


cos 7mn = 


(w m , w n ) 

Ikmll • IKII 


(73) 


The angle between any mode and itself is zero, because this formula gives 
cos7 mm = 1. The angle between modes of opposite parity is tt/2, because 
the inner product (w + m, w n ) is zero. The remaining case to consider is the 
angle between two modes of the same parity, where m+n is an even number. 
In this case, the angle is 


00S7mn 


( 1 )[ m / 2 l~H n / 2 ] COS ($ m COS^n ||VVn|| • ||^.|| 


(ll0m|| 2 + COS 2 6 m H^mll 2 ) 1/2 (||<^n|| 2 + COS 2 <5„ ||^n|| 2 ) 


2\l/2 


(74) 


The sine of the angle j mn is conveniently defined by a two-equation formula 
as follows. 


Om 


sin 7 mn 


COS Sr 


jhM 
"IhM 
1 + 0.1 


+ <£ 

, 1 + + a 2 + a 2 


.«n) 


1/2 


(75) 

(76) 


In a typical case, the norms are of the same order of magnitude, and sin 7 mn m 
y/3/2. The angle between the modes would be roughly 60°. In the case 
where k r = 1, the angles between all even modes is about 57°, and the angles 
between all odd modes is about 68°. 

The frequencies u> n associated with the cross-stream modes are 


u n = -iu(k 2 r + k 2 zn ) = -iuk 2 (77) 

Note that the frequencies are inversely proportional to the Reynolds number 
and are separated by roughly -iu(2n + 3)7r 2 /4. When the Reynolds number 
is large, there will be a considerable number of frequencies on the negative 
imaginary axis whose magnitude is small, but these frequencies are always 
distinct, that is, they form a discrete set. Waves corresponding to these 
frequencies decay slowly in time, but are absolutely stable. Note that, since 
k z0 > tt/ 2, the “closest” frequency u 0 is “below” -iu(k 2 + 7r 2 /4) on the 
imaginary axis. 
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5 Cross-Stream Spectral Method 

5.1 Standard Eigenvalue Equation 

The solutions of the cross stream propagation case may be useful in a spectral 
approximation for the general case. Since each solution satisfies the boundary 
conditions, a representation of w(z ) by an infinite sum of these terms will 
also satisfy the boundary conditions. We therefore seek an spectral solution 
to the Orr- Sommer feld equation by using the Galerkin criterion. Thus, let 

OQ 

w(z) = a n w n (z)/kn = [ w n (z)/kn J { On } (78) 

n=0 

where.it is clear that the row matrix has dimensions (l,oo) and the column 
has dimensions (oo, 1). If the sequence a n is bounded for large n, then the 
series converges as k~ 2 because the norm of w n is proportional to I fc" 1 . This 
spectral solution in the Orr-Sommerfeld equation gives an equation 

iu [ L*w n /k n J { On } + M (UL 2 - U'^Wn/kn J { 0« } = U/ L L 2 W n /k n J { fl n } 

(79) 

This equation can be simplified somewhat by using the equation for the cross 
stream modes, a special case of equation 79 which is satisfied identically for 
each term in the series. 

iu \_L Wfi/kn J { a n } = [ L 2 w n /k n J { o, n } (SO) 

Subtracting equation 80 from equation 79 gives the simplified approximation 

kg [ (UL 2 - U")w n /kn J { On } = L(w- U^tfwn/kn J { 0 „ } (81) 

Any error in this solution must be orthogonal to the elements w m (z) of the 
set of cross-stream basis functions. This gives the linear eigenvalue problem 

k x [kmiwm^-UL 2 + U")w n /k n )]{a n ) + [a>„] {a„} =u>{an} (82) 

Equation 82 can be stated concisely as 

[k x [D} + [tOn]]{an} =u{a n } (83) 
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where it is understood that [w n ] is a diagonal matrix and the matrix [D] is 
given by 

[D] = [fc m {w m ,(-^L 2 + ^> n >/^] (84) 

The above equations arc given in terms of infinite matrices, and the compu- 
tations must be restricted to finite matrices. It is assumed here, and will be 
demonstrated later, that a convergence occurs as the matrix size increases. 
But no proof of convergence will be attemped. Further discussion of the 
convergence of spectral approximations may be found in the book by F. 
Chatelin [7]. 

The matrix [ D ] is a function of k r only. Once this is constructed, the 
dependence of the eigenvalues on cos <f> can be observed by repeated solu- 
tions with different values of that parameter. The matrix [D] is real, but 
not symmetric. It is clear that the eigenvalues of equation 83, viewed as a 
function of <j>, will be stationary at the downstream direction <f> = 0. This 
is consistent with Squire’s theorem, which says that the imaginary parts of 
the eigenvalues take on extremal values at that point. The effect of Reynolds 
number is contained within the cross-mode frequencies u n . These frequencies 
are negative imaginary numbers which are inversely proportional to Reynolds 
number and directly proportional to the square of the spherical wave number. 

5.2 Evaluation of the array D mn 

The elements of the matrix will be designated by D mn , with to indicating 
the row number and n indicating the column. Using the identity 46 gives 

Dmn = km{D,Wm(^n)0Cn 4" km{U , U/ m W n ) / kn (85) 

The classic cases of Couette and Poiseuille flow are included within the case 
where U (z) is a polynomial of degree p > 2. 


p 


u(z) = J2 u j zJ 

(86) 

J=0 


P-2 


u"{z) = 5>;v 

(87) 

3=0 


u" — (j + 2 ){j + l)uj +2 

(88) 
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Equation 85 is rewritten for the polynomial velocity profile with inner prod- 
ucts of powers of z and products of the mode functions. 

P P-2 

Dmn — ^ ] Uj k m (z^ , W m (f> n) Q? n + ^ ] Uj k m {jA , W rn U} n ') / kn (89) 

J=0 j = 0 

This rewritten form shows that it is necessary only to evaluate three different 
arrays. These are 



Ajmn — {z^ j ^m^n) 

(90) 


Bjmn = 

(91) 


C jmn = (z-’.V’mV'n) 

(92) 

The other array is a transpose of the second one above 



> ifimlpn) — Bjnm — &jmn 

(93) 

Now, the array D mn is given by the long expression 


D mn = 

amain Y U 3 ( A Jmn + (-l) [m/2! COS 6 m B jmn ) 
3 = 0 


+ 

. 2 Y U 3 (Ajmn + ( l) |m/2l+[n/21 COS 6 m COS 5 n C jm „) 

j— 0 

(94) 

+ 

1.2 5Z U 'j (( 1 ) |m/2 ' 008 <5 mBjmn + ( I) 1 " 72 * COS 6 n Bj ml ? 

j = 0 


5 . 2.1 Evaluation of arrays A jTnn 

' — *. - • a •--* -- 5 • 


The elements of the array A jmn are given by the formula 



Ajmn 

(95) 


The trigonometric identities for products of trigonometric functions are used 
to convert the products to functions of the sums and differences of the re- 
spective arguments. The specific identities to be used are as follows. 


cos k zm z cos k zn z = -cos(k zm - k zn )z + -cos(k zm + k zn )z (96) 
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sin k zrn z cos h zn z ^ si n(fc 2Tn k zr P)z sin ( y k zrn -(- k zn )z (97) 

cos k ZTn z sin k zri z — sin (/c 2 m kzn)% -i- — sin(fc zm -i- k zn ^z (98) 

sin k zm z sin k zn z = ^ cos(k zm - k zn )z - ^ cos(/c zm + k zn )z (99) 

Utilizing these identities gives the following formulas for the elements A, mn . 
If m and n are both even, then 

Ajmn 2 (^ > COs(A; zm ^zn)^) 4* ^ (z J , COs(fc zm + k zn }z) (100) 
If m is odd and n is even, then 

Ajmn = -(z J ,sin(A: 2m k zn )z) + — (z^, sin(£ zm + k zn )z) (101) 
If m is even and n is odd, then 

Ajmn = ~\i z \ sin(fc 2m - k zn )z) + )^(z } , sin(fc zm + k zn )z) (102) 

Finally, if m and n are both odd, then 

Ajmn = 2 C^jCOs(fc zm k zn )z) — (z^ , cos(fc 2m + k zn ')z') (103) 


5.2.2 Evaluation of arrays B jmn 

When integrals involving ip occur, ik T would be used in place of k zm or k zn . 
This permits the development of explicit formulas for the elements Bj mn 
which are very similar to the ones for the matrix elements A jmn . If m and n 
are both even, then 


Q cos(&2n 1" 

jmU = ^shT r 


If m is odd and n is even, then 


Bjmn — 


$(z j ,sin(fc zn + iAy)z) 
sinhAv 


(104) 


(105) 
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If m is even and n is odd, then 


„ %(z j ,sm(k zn + ikr)z) 

&jmn — V I 

J cosh K r 

Finally, if m and n are both odd, then 

— cos (k zn + ikr)z) 


Bjmn — 


sinh kr 


( 106 ) 


(107) 


Note that every even-numbered row of the array Bj mn is identical, as is every 
odd-numbered row. Thus, the complete array B jmn is defined by (2 x N) 
elements. 


5.2.3 Evaluation of arrays Cjmn 

The formulas for the elements of the arrays Cj mn are found by letting both 
k xn and k zm go to ikr- If m and n are both even, then 


■'jmn 


(z j , 1) 4- (z j , cos 2 ik r z) 
1 4- cos 2 ik r 


If m is odd and n is even, then 

^ _ (z j , sin 2ikrz) 

Jmn ~ sin 2ikr 

If m is even and n is odd, then 

(z J ,sin 2ikrz) 


Cjmn — 


sin 2 ikr 


(108) 


(109) 


( 110 ) 


Finally, if m and n are both odd, then 


Cjmn 


(z J , 1) - (,2 J , COS 2ikrz) 

1 — cos 2 ikr 


(111) 


Only a (2 x 2) array must be evaluated to provide all elements of C jmn . 
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5.2.4 Evaluation of Elementary Integrals 

Let the sums and differences of the vertical wavenumbers be designated as 

&rnn (112) 

All of the inner products in the previous section are evaluated using the 
following two elementary integrals. 


{z 3 ,sma mn z) 
(z j , cos a mn z) 


( = 0, if j is even 
1 ^ 0, if j is odd 
f 7^ 0, if j is even 
1 = 0, if j is odd 


(113) 

(114) 


The parameter <r mn in equation 112 can be small or large. Diagonal elements 
in the array Aj mn have a mn = 0. The elements in the array C jmn have terms 
where <r mn = 0 or |2iAv|, both of which are small if k r is small, |fc r | << 1. 
Off-diagonal elements will generally have integrals where <r mn is large. The 
integrals where \<r mn \ is small or zero are best evaluated using power series 
for the sine and cosine functions. The following formulas presume that j is 
even or odd, depending on whether the integral has a non-zero value . 


. , “ (-]\P„ 2 P +1 

(^.sina™*) = ^ — 2*2 \ a I <• Q(\\ 

£S(2p + 2j + 2)(2p+l)!' 


(^,co s <w) = f; — 


S3(2p + 2 + l)(2p)!' 


< O(I) 


(115) 

(116) 


When |<r mn | is a finite number, then exact expressions for the inner products 
are available. 


(z * , sin a mn z) 


{ z , COS (J mn Z ) 


+ 


COS a mn ^ (~l) p j! 

&mn p=o (J 2p)!(Tmn 

Sin ^ (~1 ) P Q - 1)! 

p—0 {j 1 — 2p)!<7mri 
Sin Ornn (~l) P j! 

a mn p_o (j — 2p)!<rinn 

cos <J mn 1 (— 1) P Q - l)! 

0 -1- 2p)!<7&„ 


(117) 


(118) 
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All of the finite-series expressions above are exact, and the condition on 
|a mn | is a general recommendation for their application, but not an absolute 
requirement. The finite series are limited by machine accuracy, because their 
results are given as the differences of large numbers when |<r mn | is small. The 
two terms in the finite series are 0(|cr mn | -2 ^/ 2 ] +1 )), and their difference must 
be 0(1). Consequently, any machine error c min is amplified by the magnitude 
of the two terms. If the desired computation error is less than e, then |<j mn | 
is limited 

Cminkmnr 2(Ij/2]+1) < C (119) 

The finite series should then be used if 

K„| > (120) 

If we attempt to make the computation error as small as the machine error, 
then \<T m n\ > 1 in the finite series. Achieving machine accuracy in the infinite 
series requires enough terms in the series such that the inverse of the factorial 
is less than e mjn , that is, p ranges from 0 to M, where M is the least integer 
satisfying the condition 

(2 M + 1)! < (121) 

For example, a 64-bit computer word corresponds to a machine error of about 
10- 15 , and 19! « 1.2xl0 17 , so that 2 M + 1 = 19, or M = 9 as an upper limit 
in the infinite series with |cr mn | < 1 would be consistent with the goal of 
machine accuracy for a 64-bit word. 
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6 Application to Poiseuille Flow Stability 

This section considers the application of the proceeding analysis to the sta- 
bility of Poiseuille flow in a channel. This is the simplest case to consider 
because the matrix solutions can be made separately for the even modes and 
for the odd modes. The asymptotic behavior of the matrix elements for large 
indices gives approximate solutions for the higher modes. This approximate 
solution is given below, and then computational results will be shown for a 
test case. 

Poiseuille flow is symmetric with U{z) — 1 — z 2 . The velocity coefficients 
are then {ttj} = {1,0,— 1}. There is a single coefficient for the curvature 
of the velocity, { u '•} = {-2}. The array 4>mn is 8 mn {ot„, the array B 0mn is 
identically zero, and the array C 0m n is 6 mn Con n , so that the array D mn is 

Dmn = a m a n - A 2 mn - (-l) [m/2] cos 6 rn B 2m ^j 

-2 + (_ i)M+l»/J| cosSmCOS (^C^) (122) 

The matrices for the Poiseuille flow may be formed using only even modes or 
only odd modes, since the flow profile is symmetric. This greatly simplifies 
the computations for Poiseuille flow. All indices are then either both even 
or both odd, and functions like fl n which depend only on the parity of the 
index can conveniently be designated by a parity index p. 

p = mod (m, 2) = mod (n, 2) (123) 

6.1 Asymptotic approximations for finite kr 

Asymptotic approximations were developed previously for the vertical wave- 
numbers. Similar approximations are possible for the elements of the array 
D mn . The approximations show the properties of the general eigenvalue 
problem and will lead to an asymptotic estimate of the eigenvalues of the Orr- 
Sommerfeld equation itself. The asymptotic formulas are developed following 
a reversed order from the evaluation formulas in the previous section, that is, 
the formulas for the elementary functions are developed first, then formulas 
for elementary integrals are given, and finally, the formulas for the arrays are 
enumerated. 
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6.1.1 Asymptotic formulas for wavenumber variables 

All formulas to be given here assume that K- is finite. We use less accurate 
formulas than the ones given previously for the vertical wavenumbers because 
the point of this development is to show general properties rather than to 
develop an accurate computational procedure. 

The offset variable 6 n is given by a one- term approximation 

6n~ — + , K n »kr (124) 

/Cn 

The vertical wavenumbers are given by two-term approximations. 


k zn ~ K-n 



Kn» kr 


(125) 


The spherical wavenumbers kr » depend on the radial wavenumber and the 
vertical wavenumbers. 


kn K n 



Kn» kr 


(126) 


This result gives a simple formula for the characteristic frequencies of the 
cross-modes 


UJ n 


~ —iVK 


2 

n 



/Cf| >> kr 


(127) 


The mode normalization constants are also easy to estimate with these re- 
sults. 


an ~^[ l + £iJ + 


Kn » K 


(128) 


6.1.2 Elementary integrals 

The elementary integrals with small arguments can occur only when the 
terms represent positions on the diagonals of the arrays and when k± is the 
difference of the wavenumbers. But in this case, the arguments are exactly 
zero. The elementary integrals are then 

{z , cos(fc zm k zn )z) = TYi = n (129) 
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In the more general case, the arguments are finite. 

2 &mn ^ 0 


/ 2 \ ^n (T mn ( 2 \ COS <T mn 

(z, COS <T mn z) = 1 — + 2 - 

<^mn \ ff mn / 


(130) 


The array Comn is given by 


1 


a_- 


(131) 


6.1.3 Asymptotic formulas for arrays 

The asymptotic form for the array B 2mn shows that it is inversely propor- 
tional to k 2 . 

B 2 mn ~ A(-l)[("+2)/ 2 )] + ... ( 132 ) 

K n 

The general form of the array D mn can be completed with as estimate of 
the array A 2mn . Two cases must be considered here. In the first case, the 
indices are equal, and the elements are of order unity. In the other case, the 
indices are unequal, and the elements are small. The case of equal scripts is 
given by 

A * nn ~ 6 + ( _1 ) n 4/c2^ + '“ ( 133 ) 

When the scripts are unequal, we use the exact equation 113 with the asymp- 
totic form for the sum and difference of the offsets. 



If m and n are nearly equal, that is, the elements are near the principal 
diagonal of the array, then this approximation is of order («)~ 2 when the 
difference is taken and of order (/c) _I when the sum is taken. In either case, 
the next term is of order (/c) -3 , so that it consistent to use it in both cases 
with this understanding. This gives the following result for the difference of 
the wavenumbers. 


kzm i k zn ^ (ft m i AC n ) ( 1 ^ 


(■ 



«m«n 



(135) 
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These formulas give an asymptotic form for the inner product of z 2 and the 
cosine functions. First, express the exact integral in terms of the wavenum- 
bers and offset variables. 


{z 2 ,cos(k zm ± k^z) = -(-1) 


win sin(6 m ± 6 n ) 


+ 2 (— 1 ) 


(k zm i kzn) 

min COS (6m i 6n) 


(k zm ± k zn ) 2 

The sine and cosine of the offset variable pairs are given by 


sin(6 m ±6„) ~ ±(/c m ±/Cn) 


(— ) 


+ 


/c c x , (K m ±K n ) 2 ( 0 p \ , 

~ 1 J + 

These formulas are used to give the terms in equation 136. 
sin(6 m ± <5„) 


l-) + - ' 

__j (Jk\ 

{k Z m ^ k zn ) 3 (K m ± K n ) 2 \K m K n J 


(kzm i kzn) 
sin(6 m ± 8 n ) 


± 

±- 


+ 


cos(6 m ± 6n ) 

(k zm ± k zn ) 2 ( 


K m i K n ) \ K,m^-n J 


+ 


(136) 


(137) 

(138) 

(139) 

(140) 

(141) 


The final asymptotic formula for the inner products of z 2 and the cosine 
functions is 


(z 2 ,COs(k Z m±kzn)z) ~ (- 1 )^" 7 ~ | x 2 T 

This gives a simple result for the dominant term in the array A^mn- 

1 


(142) 


*2mn 


(-I)' 


(K m K n ) 2 


(143) 


32 



6.1.4 Asymptotic formula for array D mn 

The above results produce an easy formula for the array D mn which shows 
that the matrix is dominated by the diagonal elements. That is, the off- 
diagonal elements in any given row or column are smaller in magnitude than 
the diagonal element. The matrix is clearly dominated by the first line in 
equation 94 because the terms in the second line vary as k~ 2 . In the first line, 
terms from the array B 2mn also vary as k ~ 2 so that they will be negligible 
in comparison to the dominant terms in A 2 mn , which are of order one near 
the diagonal. Thus, a first-order estimate of the array D mn for large row and 
column indices m, n is 


D L 


■( 


(-D 1 


2 


, if n = m; 
, if n ^ m. 


(144) 


Since m and n have the same parity, they must differ by some multiple of 
2. The elements of the diagonals adjacent to the principal diagonal then 
have magnitude 2/7r 2 , or about 1/5. The diagonal elements are all nearly 
2/3. This is the average value of the flow speed U(z) for the present case of 
Poiseuille flow. The symbol T mn is used because the matrix elements depend 
only on position with respect to the principle diagonal. Matrices of this type 
are called Toeplitz matrices. The matrix D mn is asymptotic, when both row 
and column indices are large, to the Toeplitz matrix T mn . 


6.2 Gersgorin region and asymptotic eigenvalues 

It is known from matrix theory [9] that the eigenvalues of a matrix must 
lie within a domain in the complex plane called the Gersgorin region. This 
region is the union of all points which lie within circular subdomains called 
Gerggorin disks. The center of each disk is a point given by a di ag onal 
element of the matrix, and the radius of each disk is the sum of the absolute 
values of the off-diagonal elements in the corresponding row or column. This 
actually defines two Gersgorin regions, but, in the present case of a Toeplitz 
matrix, the row and column sums are identical so that the row-disks and 
column-disks are identical. The Toeplitz matrix here is particularly simple 
and the disk radii are bounded by an infinite series with known sum. This 
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bound is given by 


K = 


2| k x 

7T 2 


d? 0 


d 2 


< 


4|fc* 


7T* 


°° 1 2 


(145) 


That is, each GerSgorin disk has the same radius, and this radius is less than 
the product of the average flow speed and the axial wavenumber. In the 
limiting case of cross-modes, the axial wavenumber approaches zero and the 
eigenvalue must approach the disk center, or the value of a diagonal element 
of the matrix. 

The center of the Gersgorin disk is a good first guess for the eigenvalue 
even when the axial wavenumber is not zero, at least in the case of large 
mode index n where the disks are disjoint regions. 


Wn « = k x D nn - ivkl ~ -k x 


U'KZ 


(146) 


The real part of the approximate eigenvalue (disk center) is the product of 
the average flow speed and the axial wavenumber k x = kj. cos<f), and the 
imaginary part of the wavenumber is identical to the complex frequency of 
the cross-stream mode. 

Figure 9 depicts two Gersgorin disks as defined by the above equations. 
An eigenvalue, called u n must lie within the shaded domain centered at £>„. 
Since the disk radius is less than to the real part of the diagonal element, 
there is no possibility that the real part of the eigenvalue is negative when 
the axial wavenumber is positive. More generally, the sign of the real part of 
the eigenvalue must match the sign of the axial wavenumber k x . 



if k x > 0; 
if k x < 0. 


(147) 


The imaginary part of the eigenvalue may be positive or negative if the index 
n, is small, but, if 

vnl > (148) 

then the imaginary part of the eigenvalue would be required to be negative, 
and all eigenvalues meeting this criterion must then lie in the lower half of 
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the complex plane and represent stable modes. Modes with index n are then 
stable if 

n > 2 

The Gersgorin disks may intersect as shown in figure 9, but when the mode 
index is large, the disks become disjoint and each disk must contain a single 
eigenvalue. The condition for disjoint disks is 



n > 


/4 Re 
V 37r 2 


3 ) 


(150) 


The Gersgorin criterion guarantees only that an eigenvalue lies within the 



Figure 9: Gersgorin disks and eigenvalues for Poiseuille flow. The large 
circles are the outer bounds for the disks. Shaded areas are actual disks. 
Disk centers are open circular symbols and eigenvalues are filled circular 
symbols. 


disk, not where it lies, but the computations in the following section will 
demonstrate that the eigenvalue approaches the disk center when the index is 
large. That is, the “higher” eigenvalues approach the centers of the GersSgorin 
disks. 
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6.3 Asymptotic convergence of eigenvalues 

A simple approximation to the general eigenvalue problem for Poiseuille flow 
is given by the asymptotic eigenvalue problem. 

[k x [T mn \ - iv [«*] ] { a n } = u { an } (151) 

The matrices in this equation are within a few percent of the actual matrices 
so that this equation contains the principle features of the actual problem. 
This similarity is used here to demonstrate the convergence of the eigenvalues 
to the asymptotic formula. This is done by comparisons of results from the 
asymptotic formula and computations from 32nd order, 48th order, and 64th 
order matrices. 



Figure 10: Even eigenvalues iu for Poiseuille flow based on a 32nd order 
matrix, Av = 1,^ = 0, and Re = 10 4 . Solid symbols represent computed 
eigenvalues using the Toeplitz matrix [T] and the asymptotic values of the 
cross-mode wavenumbers, K n . Gray symbols represent the asymptotic for- 
mula for the eigenvalues. 

Results from the 32nd order matrix solution are shown in figure 10. The 
32 complex eigenvalues are distinct, and fall within a unit square in the 4th 
quadrant of the complex plane. This is a stable quadrant, according to the 
sign convention adopted for this paper. This approximate computation indi- 
cates flow stability, whereas it is known that this flow is theoretically unstable 
for Reynolds numbers above about 5772. But the purpose of these results 
is to show comparisons to the asymptotic formula, not make an absolutely 
accurate computation. The eigenvalues predicted by the asymptotic formula 
lie on a vertical line in the complex w-plane which is defined by the average 
flow speed, 2/3. Since the length of this line gets large as more eigenvalues 
are added, the plot shows iu to make the figure fit more comfortably on the 


36 


page. Note that only a few of the computed eigenvalues fall near the points 
generated by the asymptotic formlula. This is an indication that the selected 
matrix order is not sufficiently large. 



Figure 11: Even eigenvalues iu for Poiseuille flow based on a 48th order 
matrix, kr = 1, <f> — 0, and Re = 10 4 . Symbols are the same as in figure 10. 

The results of computations using a 48th order matrix are shown in figure 
11. It is seen that the more-attenuated eigenvalues have moved in the direc- 
tion of the asymptotic formula results, while the less-attenuated eigenvalues 
are essentially unchanged. As expected, the increase in matrix size appears 
to produce a greater number of more-accurate eigenvalues. Eigenvalues from 



Figure 12: Even eigenvalues iu for Poiseuille flow based on a 64th order 
matrix, Av = 1, <f> — 0, and Re = 10 4 . Symbols are the same as in figure 10. 

the 64th order matrix are shown in figure 12. It is seen that all but about 
6 have approached the asymptotic results. It appears that the computed 
spectrum is approaching the asymptotic formula in the domain where the 
mode index n is large. 
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6.4 Computational results for Poiseuille flow. 

The previous results show the qualities of the computation using an asymp- 
totic approximation to the matrix eigenvalue problem. This section will show 
results of the actual matrix eigenvalue problem. 

6.4.1 Numerical conditioning 

The matrix eigenvalue problem with the cross-stream basis is well-condition- 
ed in comparison to the Chebyshev/Tau method. This is demonstrated in 
figure 13 for the same physical case as in figure 3. That is, computations 
were made for the even modes using an order 32 matrix. As in figure 3, 
the computations were made with 32-bit and 64-bit words as a check of 
the sensitivity. Most of the 64-bit (black) eigenvalues plot directly onto the 



Figure 13: Even eigenvalues for Poiseuille flow using the Cross-Stream 
method and order 32 matrices. Symbols and physical parameters are the 
same as in figure 3. 

the 32-bit eigenvalues, demonstrating the relative insensitivity to computer 
word size. More importantly, a careful count of the symbols shows that all 
32 eigenvalues are accounted for within the displayed portion of the complex 
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plane. The eigenvalues are distinct, and a sorting operation on the imaginary 
part can then select a “least-attenuated” mode without the ambiguity of a 
partial eigenvalue list. There is, however, a discemable effect of word size on 
a few of the eigenvalues, so that the remainder of the computed results will 
be based on 64-bit words. 


6.4.2 Asymptotic convergence of eigenvalues 


Figure 14 shows the computations for a 64th order matrix. Again, gray 
symbols are the centers of the Gersgorin disks (D nn - ti/fc*) whereas the 
computed results are shown as black symbols. The radii of the Gersgorin 
disks will be roughly (2/3 )k x so that the higher eigenvalues probably lie in the 
4th quadrant. Inspection of figure 14 indicates that the computed eigenvalues 
lie even closer to the Gersgorin centers than in the case of the approximate 
matrix eigenvalue problem. Camparing figures 12 and 14 indicates that the 
Toeplitz matrix approximation affects primarily the lesser-attenuated modes 
which propagate at less than the average flow speed. 
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Figure 14: Even eigenvalues for Poiseuille flow using the cross-stream basis 
and order 64 matrices. Symbols and physical parameters are the same as in 
figure 12. 


Figures 15, 16, and 17 show the convergence of the eigenvalues with increasing 
matrix size. The matrix order is doubled in each successive figure, starting 
from figure 15 where the order is 64, so that figure 16 has matrix order 128 
and figure 17 has matrix order 256. figure 15 is the same 64-bit data as in 
figure 14, but only distance of the eigenvalue from the Gersgorin center is 
plotted. The abscissa in each figure is the actual mode index, which ranges 
up to 510 in figure 17. These computations show clearly that the eigenvalues 
converge to the Gersgorin centers as the index increases. It can be concluded 
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from this that the cross-stream method provides a means to compute an 
indefinite number of eigenvalues of the Orr-Sommerfeld equation. The first 
several hundred are computable by standard matrix eigenvalue algorithms, 
and any countable number of higher eigenvalues are approximated by the 
easily-computable Gersgorin centers. 
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Figure 15: Distance of even eigenvalues from Gerggorin centers based on 
matrix order 64 computations. Physical parameters are the same as in figures 


10 - 12 . 
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Figure 16: Distance of even eigenvalues from Gersgorin centers based on 
matrix order 128 computations. Symbols and physical parameters are the 
same as in figure 15. 



Figure 17: Distance of even eigenvalues from Gersgorin centers based on 
matrix order 256 computations. Symbols and physical parameters are the 
same as in figure 15. 
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7 Concluding Remarks 

The Chebyshev polynomial basis for spectral solutions to the Orr-Sommerfeld 
equation leads to a generalized matrix eigenvalue problem. Matrices in this 
problem are ill-conditioned so that solutions require large computer words. 
The conditioning problem is exacerbated with increasing matrix order so that 
the Chebyshev basis can produce only a limited number of meaningful eigen- 
values. The Chebyshev basis also produces eigenvalues in all four quadrants 
of the complex plane, even for cases believed to be stable, so that a definition 
of “least attenuated” mode is ambiguous. 

Exact solutions to the Orr-Sommerfeld equation are available for waves 
propagating perpendicular to the flow direction. These solutions, called 
cross-stream modes here, are a countable set of linearly independent func- 
tions which serve as a basis for a Hilbert space. The cross-stream modes 
are not orthogonal, but may be subdivided into sets of even and odd func- 
tions which are, or course, orthogonal. The “angles” between each of the 
even functions are roughly 57°, and the angles between the odd functions 
are about 68°. The Hilbert space with this basis serves well to estimate the 
spectrum of the Orr-Sommerfeld operator and its boundary conditions. 

The cross-stream basis leads to a standard matrix eigenvalue problem, 
which is simpler to solve than the generalized problem produced by the 
Chebyshev basis. Conditioning of this matrix eigenvalue problem is signifi- 
cantly better than the Chebyshev problem. The matrix is strictly diagonally 
dominant, and all matrix elements, for polynomial flow profiles, are given by 
known elementary integrals. 

The matrix in the cross-stream eigenvalue problem approaches the sum 
of a constant real Toeplitz matrix and a complex diagonal matrix when the 
mode index is large. The complex diagonal represents the purely damped 
frequencies of the cross-stream modes, that is, the elements of this diagonal 
matrix are negative imaginary numbers. 

The eigenvalues of the cross-stream matrix lie within well-defined Gers- 
gorin disks. It has been shown that the radius of each disk, is less than about 
2k x /3 for Poiseuille flow, where the fraction 2/3 represents the average speed 
of the flow and k x is the given axial wavenumber. As the wavenumber de- 
creases to zero, the disks become points and the eigenvalues are exactly the 
eigenvalues of the cross-stream modes. The Gersgorin disks for the cross- 
stream matrix show that the real part of all eigenvalues must have the same 
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algegraic sign as the given axial wavenumber. The Gerggorin disks show also 
that only a finite number of modes may be unstable, that is, have eigenval- 
ues with a positive imaginary part. This number increases in proportion to 
the square root of the Reynolds number. An asymptotic formula has been 
developed for the “higher” eigenvalues, which are those with greater atten- 
uation. The real part of the asymptotic eigenvalue is the product of the 
average flow speed and the axial wavenumber, and the imaginary parts are 
the negative-imaginary frequencies of the cross-stream modes. 

Computations with the cross-stream method have been made for the well- 
studied case of a Poiseuille flow with Reynolds number 10,000. Matrix orders 
up to 256, producing 512 eigenvalues, half even and half odd, show that the 
lower eigenvalues match the results of the Chebyshev method to eight signif- 
icant figures while the higher eigenvalues approach the asymptotic formula 
to within three significant figures. 

Computations also show only a single eigenvalue with positive imaginary 
part. This is the one corresponding to the critical mode defined by other 
investigators. All other eigenvalues have negative imaginary parts. The 
cross-stream method permits an unambiguous definition of “least attenuated 
/ mode” through a sorting process on the imaginary part on all eigenvalues 

produced by the computation. 

The cross-stream method is a well-conditioned and robust computational 
approach which can produce an essentially unlimited number of accurate 
eigenvalues of the Orr-Sommerfeld equation. 
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